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Abstract. This paper concerns with numerical approximations of solutions of second order 
fully nonlinear partial differential equations (PDEs). A new notion of weak solutions, called moment 
solutions, is introduced for second order fully nonlinear PDEs. Unlike viscosity solutions, moment 
solutions are defined by a constructive method, called vanishing moment method, hence, they can 
be readily computed by existing numerical methods such as finite difference, finite element, spectral 
Galerkin, and discontinuous Galerkin methods with "guaranteed" convergence. The main idea of the 
proposed vanishing moment method is to approximate a second order fully nonlinear PDE by a higher 
order, in particular, a fourth order quasilinear PDE. We show by various numerical experiments 
the viability of the proposed vanishing moment method. All our numerical experiments show the 
convergence of the vanishing moment method, and they also show that moment solutions coincide 
with viscosity solutions whenever the latter exist. 
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1. Introduction. Fully nonlinear PDEs are those PDEs which depend nonlin- 
ear ly on the highest order derivatives of unknown functions. Fully nonlinear PDEs 
arise from many areas in science and engineering such as kinetic theory, materials sci- 
ence, differential geometry, general relativity, optimal control, mass transportation, 
image processing and computer vision, meteorology, semigeostrophic fluid dynamics. 
They constitute the most difficult class of differential equations to analyze analytically 
and to approximate numerically, see [T3| [46l [45l [151 EHl El] and references therein. 

The general first order fully nonlinear PDE has the form 

F{\/u{x), u{x),x) = X encIC. (1.1) 

The best known examples include Eikonal equation 

\Vu{x)\ = f{x) xeQ, 

and the general Hamilton- Jacobi equation [46l [25] 

H{\/u{x)) = xen. 

The general second order fully nonlinear PDE, which will be the focus of this paper, 
takes the form 

F{D^u{x), \/u{x), u{x),x) = X en. (1.2) 
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where and throughout this paper D'^u{x) denotes the Hessian matrix of u at x. The 
best known examples are the Monge- Ampere type equations |46j |511 157] 

det{D'^u{x)) = f{\/u{x), u{x),x) x e 

and the Behman equations |46l [45] 

sup Lq{D'^u^ Vii, x) = 0, (1-3) 

eee 

where det{D'^u{x)) stands for the determinant of D'^u{x) and Lq is a given family of 
second order linear differential operators. 

For the first order fully nonlinear PDEs, tremendous progresses have been made in 
the past three decades. A revolutionary viscosity solution theory has been established 
(cf. [26l |24l |25l |45]) and wealthy amount of efficient and robust numerical methods 
and algorithms have been developed and implemented (cf. [TOl l2Ql l27l l58l l67l ITTl [72] ) . 
However, for second order fully nonlinear PDEs, the situation is strikingly different. 
On one hand, there have been enormous advances in theoretical analysis in the past 
two decades after the introduction of the notion of viscosity solutions by M. Crandall 
and P. L. Lions in 1983 (cf. [H [H [H [25l|46l [51] ) • On the other hand, in contrast 
to the success of the PDE analysis, numerical solutions for general second order fully 
nonlinear PDEs (except in the case of Bellman type PDEs, see below for details) is 
mostly an untouched area, and computing viscosity solutions of second order fully 
nonlinear PDEs has been impracticable. There are several reasons for this lack of 
progress. Firstly, the strong nonlinearity is an obvious one. Secondly, the conditional 
uniqueness (i.e., uniqueness holds only in certain class of functions) of solutions is 
difficult to handle numerically. Lastly and most importantly, the notion of viscosity 
solutions, which is not variational, has no equivalence at the discrete level. 

To see the above points, let us consider the following model Dirichlet problem for 
the Monge- Ampere equation: 

det{D^u) = f inn, (1.4) 
u = g on dft. (1.5) 

It is well-known that for non-strictly convex domain Q the above problem does not 
have classical solutions in general even /, g and dft are smooth (see [46 ). Classical 
result of A. D. Aleksandrov states that the Dirichlet problem with / > has a 
unique generalized solution in the class of convex functions (cf. [TJ [17]). Major 
progress on analysis of problem ( |1.4| )-( |1.5[ ) has been made later by using the viscosity 
solution concept and machinery (cf. [13( I25( [SP). We recall that a convex function 
u G C^{n) satisfying u = g on dft is called a viscosity subsolution (resp. viscosity 
supers olution) of ( |1.4[ ) if for any cp e C'^ there holds det{D'^(p{xo)) < f{xo) (resp. 
det{D'^(f{xo)) > f{xo)) provided that u — (f has a local maximum (resp. a local 
minimum) ai xq G Q. u G C^{Q) is called a viscosity solution if it is both a viscosity 
subsolution and a viscosity supers olution. 

First, the reason to restrict the admissible set to be the set of convex functions 
is that the Monge- Ampere equation is elliptic only in that set [51 . It should be 



noted that in general the Dirichlet problem (1.4)-(1.5) may have other (nonconvex) 
solutions besides the unique convex solution, multiple solutions are often expected for 
the Monge- Ampere type PDEs and for second order fully nonlinear PDEs. It is easy 



to see that if one discretizes (1.4) straightforwardly using the finite difference method. 
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one immediately loses control on which solution the numerical scheme approximates 
even assuming that the nonlinear discrete problem has solutions. Second^ the situation 
is even worse if one tries to formulate a Galerkin type method (such as the finite 
element method and the spectral Galerkin method), because there is no variational 
or weak formulation to start with. In fact, this is clear from the definition of viscosity 
solutions. It is not defined by the traditional integration by parts approach, instead, it 
is defined by a ^^differentiation by parts^^ (a terminology coined by L. C. Evans [33 , 34]) 
approach. Although the ^^differentiation by parts^^ approach has worked remarkably 
well for establishing the viscosity solution theory for second order fully nonlinear 
PDEs in the past two decades, it is extremely difficult (if all possible) to mimic it 
at the discrete level. Third, regardless which method is used, one can easily envisage 
that the anticipated algebraic problem from the discretization of a fully nonlinear 
PDE such as the Monge- Ampere equation must be very difficult to solve due to the 
nonuniqueness of solutions and very strong nonlinearity. 

Nevertheless, a few recent numerical attempts and results have been known in 
the literature. In [66] Oliker and Prussner proposed a finite difference scheme for 
computing Aleksandrov measure induced by D'^u (and obtained the solution u of 
(1.4) as a by-product) in 2-d. The scheme is extremely geometric and difficult to 
use and to generalize to other second order fully nonlinear PDEs. In [7 Barles and 
Souganidis showed that any monotone, stable and consistent finite difference scheme 
converges to the correct solution provided that there exists a comparison principle for 
the limiting equation. Their result provides a guideline for constructing convergent 
finite difference methods although it did not address how to construct such a scheme. 
Very recently, Oberman [65 was able to construct some wide stencil finite difference 
schemes which fullfil the criterions listed in [7 for the Monge- Ampere type equations. 
In [5] Baginski and Whitaker proposed a finite difference scheme for Gauss curvature 
equation (see Q in 2-d by mimicking the unique continuation method (used to prove 
existence of the PDE) at the discrete level. Finally, in a series of papers [28^1291 130 1 131 ] 
Dean and Glowinski proposed an augmented Lagrange multiplier method and a least 
squares method for problem (1.4)-(1.5) and the Pucci's equation (cf. [TT, 'W) in 2-d 
by treating the Monge-Ampere equation and Pucci's equation as a constraint and 
using a variational criterion to select a particular solution. Numerical experiments 
results were reported in [66^, |65l [S] [28l |29l |30l |3T] , however, convergence analysis was 
not addressed except in [65 . 

In addition, we like to remark that there is a considerable amount of literature 
available on using finite difference methods to approximate viscosity solutions of sec- 
ond order fully nonlinear Bellman type PDE (1.3) arisen from stochastic optimal 
control. See [Zl El iMl IS^ • Due to the special nonlinearity of the Bellman type PDEs, 
the approach used and the methods proposed in those papers unfortunately could not 
be extended to other types of second order fully nonlinear PDEs since the construction 
of those methods critically relies on the linearity of the operators Lq. 

The first goal of this paper is to introduce a new weak solution concept and 
a method to construct such a solution for second order fully nonlinear PDEs, in 
particular, for the Monge-Ampere type equations. These new weak solutions are 
called moment solutions and the method to construct such a moment solution is called 
the vanishing moment method. The crux of this new method is that we approximate 
a second order fully nonlinear PDE by a sequence of higher order (in particular, 
fourth order) quasilinear PDEs. The limit of the solution sequence of the higher 
order PDEs, if exists, is defined as a moment solution of the original second order 
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fully nonlinear PDE. Hence, moment solutions are constructive by nature. The second 
goal of this paper is to present a number of numerical methods for computing moment 
solutions of second order fully nonlinear PDEs, and to present extensive numerical 
experiment results to demonstrate the convergence and effectiveness of the proposed 
vanishing moment methodology. Indeed, one of advantages of the vanishing moment 
method is that it allows one to use wealthy amount of existing numerical methods 
and algorithms as well as computer codes for fourth order linear and quasilinear 
PDEs to solve second order fully nonlinear PDEs. The third and the last goal of 
this paper is to show using numerical studies that the notion of moment solutions 
generalizes the notion of viscosity solutions in the sense that the former coincides 
with the later whenever the later exists. These numerical studies indeed motivate 
us to give a rigorous convergence analysis of the vanishing moment method for the 
Monge- Ampere equation in two spatial dimensions [36] . 

The remainder of the paper is organized as follows. In ^ we introduce the 
abstract framework of moment solutions and the vanishing moment method for gen- 
eral second order fully nonlinear PDEs. In ^ we propose two classes of numerical 
discretization methods and briefly discuss solution algorithms. In Q we apply the ab- 
stract framework to several classes of second order fully nonlinear PDEs which include 
the Monge-Ampere type equations, Pucci's extremal equations, the infinite Laplace 
equation, and second order parabolic fully nonlinear PDEs. In ^ we present many 
2-d and 3-d numerical experiment results to demonstrate the convergence and effec- 
tiveness of the vanishing moment methodology, and provide numerical evidences of 
the agreement of moment solutions and viscosity solutions whenever the latter exists. 
The paper is concluded by a summary and some conclusions in §6. 

2. Vanishing moment method and the notion of moment solutions. 

2.1. Prehminaries. Standard space notation will be adopted throughout this 
paper, we refer to [46l [57] for their exact definitions. Q denotes a generic bounded 
domain in R^. (•, •) and (•, •) are used to denote the L^-inner products on Q and on 
dQ, respectively. We assume n > 2, except in ^ and ^ where we restrict n = 2, 3 
when we develop numerical methods and perform numerical experiments. 

Since the notion of viscosity solutions has been and will continue to be referred 
many times, and it is closely related to the notion of moment solutions to be described 
later in this paper, for readers' convenience, we briefly recall its definition and history 
here and refer to J3l |25l [26l |34] for detailed discussions. 

Definition 2.1. Suppose F : R"' x R x R"' H is continuous (nonlinear) 
function. 

(i) A function u G C^{yt) is called a viscosity subsolution of if, for every 

function ip = (fix) such that u — cp has a local maximum at x^ G ft, there 
holds 

F{\/f{x^),f{x^),x^)<0. 

(a) A function u G C^{Q) is called a viscosity supersolution of if, for every 
function (f = (p{x) such that u — cp has a local minimum at x^ G ft, there 
holds 



(Hi) A function u G C^{Ct) is called a viscosity solution of (1.1) if it is both a 
viscosity subsolution and a viscosity supersolution. 
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It should be pointed out that the above definition is a modern definition of vis- 
cosity solutions for It can be regarded as a "differentiation by parts" definition 
(cf. [34 ). However, viscosity solutions were first introduced differently by a vanishing 
viscosity procedure (cf. [26 ), that is, equation is approximated by the second 
order quasilinear PDEs 

-eAu' + F{Vu', u\ x) = 0, 

and lim^^i^, if exists, is called a viscosity solution of It was later proved that 

the two definitions are equivalent for equation (cf. [24]). 

Another important reason to favor the modern "differentiation by parts" defini- 
tion is that the definition and the notion of viscosity solutions can be readily extended 
to second order fully nonlinear PDEs. 

Definition 2.2. Suppose F : R^'^^^xR^xRxR"' ^His continuous (nonlinear) 
function. 

(i) A function u G C^{yt) is called a viscosity subsolution of \1.2\ if, for every 
function ^ = (p{x) such that u — (p has a local maximum at x^ G there 
holds 

F(I)V(^°), V(/:)(x^),(/:)(xO),xO) < 0. 



(a) A function u G C'^iQ) is called a viscosity supersolution of (1.2) if, for every 
function if> = (f>{x) such that u — cp has a local minimum at x^ G ft, there 
holds 

F(I)V(^°),V(/p(xO),(^(xO),xO) >0. 



(Hi) A function u G C^{Q) is called a viscosity solution of (1.2) if it is both a 

viscosity subsolution and a viscosity supersolution. 
As it is known now, a successful theory of viscosity solutions has been established 
for second order fully nonlinear PDEs in the past two decades (cf. [131 [ISl [51]). 
On the other hand, it should be noted that the phrase "viscosity solution" loses its 
original meaning in this theory since it has nothing to do with the vanishing viscosity 
method in the case of second order fully nonlinear PDEs. We recall that to establish 
the existence of viscosity solutions the technique used to substitute for the vanishing 
viscosity method in the theory is the classical Perron's method (cf. [25l[T3]). To the 
best of our knowledge, viscosity solutions of second order fully nonlinear PDEs were 
never defined and/or constructed by a limiting process like one described above for 
the Hamilton-Jacobi equation. 

2.2. General framework of the vanishing moment method. For the rea- 
sons and difficulties explained in Q as far as we can see, it is unlikely (at least 
very difficult if all possible) that one can directly approximate viscosity solutions of 
general second order fully nonlinear PDEs such as Monge- Ampere type equations 
using any available numerical methodology (finite difference method, finite element 
method, spectral method, meshless method etc.). From computational point of view, 
the notion of viscosity solutions is a "bad" notion for second order fully nonlinear 
PDEs because it is not constructive nor variational, so one has no handle on how to 
compute such a solution. 

In searching for a "better" notion of weak solutions for second order fully nonlinear 
PDEs, we are inspired by the following simple but crucial observation: the essence 
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of the vanishing viscosity method for the Hamilton- Jacobi equation and the original 
notion of viscosity solutions is to approximate a lower order fully nonlinear PDE by 
a sequence of higher order quasilinear PDEs. This observation then suggests us to 



apply the above principle to second order fully nonlinear PDE (1.2), this is exactly 
what we are going to do in this paper. That is, we approximate equation (1.2) by the 
following higher order quasilinear PDEs: 

GsiD'^u') + F{D^u', V^x^ x) = (r > 3, £ > 0), (2.1) 

where {G^} is a family of suitably chosen linear or quasilinear differential operators 
of order r. The above approximation then naturally leads to the next definition. 



Definition 2.3. Suppose that solves (2.1) for each e > 0, we call lim 

^ ^ £^0 + 



a moment solution of (1.2) provided that the limit exists. We also call this limiting 
process the vanishing moment method. 

Clearly, the above definition is a loose definition since the operator is not 
specified, nor is the meaning of the limit, but they will become clear later in this 
section. We note that the reason to use the terminology "moment solution" will 
also be explained later in this section, and the notion of moment solutions and the 
vanishing moment method are clearly in the spirit of the (original) notion of viscosity 
solution and the vanishing viscosity method [26 . 

To establish a complete theory of moment solutions and vanishing moment method 
for second order fully nonlinear PDEs, there are many issues we must address. For 
instance, 

• How to choose the operator Gs ? 

• What additional boundary condition(s) should satisfy? 

• Does the limit lim always exist? If it does, what is the rate of convergence? 

• How do moment solutions relate to viscosity solutions? 

• How to solve (2.1) numerically? 



• Error estimates, nonlinear solvers, computer implementations. 
As expected, we do not have answers for all the questions now, nor do we intend to 
address all of them in this paper. Instead, the focus of this paper is to develop the 
framework for moment solutions and the vanishing viscosity method, and to present 
numerical evidences to show effectiveness of the method and to justify the proposed 
approach. On the other hand, we do plan to address all theoretical issues in forth- 
coming papers |36l ES] . 

Regarding to the first issue, although the choices for G^ are abundant and flexible, 
the following are some guidelines for choosing a good operator G^- 

(a) Gs must be a linear or quasilinear operator. 

(b) ^ in some reasonable sense as e ^ 0+. 

(c) Gs{D^u) is better to be elliptic, in particular, when PDE ( |1.2[ ) is elliptic. 



(d) Equation (2.1) should be relatively easy to solve numerically. 



Since an elliptic operator is necessarily of even order, so guideline (c) above implies 



that r must be an even number in (2.1). Hence, the lowest order of equation (2.1) is 
r = 4. When talking about fourth order elliptic operators, the biharmonic operator 
stands out immediately. So we let 

Gs{D^v) := -eA^v, 



then equation (2.1) becomes 

- sA^u' + F{D'^u', V^x^ x) = 0. (2.2) 
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After the differential operator is chosen, next we need to take care the boundary 



conditions. Here we only consider Dirichlet problem for (1.2). Suppose that 



u ■■ 



on 



it is obvious that we need to impose 

— g or ^ g on dCl. 



(2.3) 



(2.4) 



Moreover, since (2.2 ) is a fourth order PDE, in order to uniquely determine we need 
to impose an additional boundary condition for u^. Mathematically, many boundary 
conditions can be used for this purpose. Physically, any additional boundary condition 
will introduce a "boundary layer" , so a better choice would be one which minimizes the 
boundary layer. Here we proposed to use one of following three boundary conditions 



or 



dAu' 
dn 



or D^u^n n = 



on dQ. 



(2.5) 



In particular, the first two boundary conditions, which are natural boundary condi- 
tions, have an advantage in PDE convergence analysis [36l[63]. Another valid bound- 
ary condition is the Neumann boundary condition = 6^ on dQ. But since this is 
an essential boundary condition, it produces a larger boundary layer than the above 
three boundary conditions. The rationale for picking the above boundary conditions 
is that we implicitly impose an extra boundary condition Au^ + li^ = ^ -|- 5^+2 on 
dVt for equation (2.2), which is a higher order perturbation of the original Dirichlet 
boundary condition = g ov^ dQ. Intuitively, we expect and hope that the ex- 
tra boundary condition converges to the original Dirichlet boundary condition as e 
tends to zero. We note that m can be any positive integer, and power 2 is used for 
convenience and it can be replaced by any positive integer. 

We now remark that when n = 2 in mechanical applications often stands for 
the vertical displacement of a plate and D'^u^ is the moment tensor^ and in the weak 
formulation, the biharmonic term becomes —e{D'^u^^D'^v) which should vanish as 
£ ^ 0+. This is the very reason why we call lim ^ if exists, a moment solution and 

call the limiting process the vanishing moment method. 

In summary, we propose to approximate the second order fully nonlinear Dirichlet 



problem (1.2), (2.3) by the fourth order quasilinear boundary value problems (2. 2), (2. 4), 
(2.5). Since we expect G W^'^ip) for m > 2,p > 2, so the convergence lim 

can be understood in i^^-topology or in i^^-topology or even in 



2.3 



in Definition 

L^-topology. To distinguish these different limits, we introduce the following refined 
definition of Definition 12.31 

Definition 2.4. Suppose that e H'^{Q) solves problem ([2^, ([2^, ([2^). 



lim is called respectively a sub- weak, weak and strong moment solution to problem 

£ ->0 + 

(1.2)^(2.3) if the convergence holds in L'^-, H^- and H'^ -topology. 

Remark 2.1. Since sub-weak and weak moment solutions do not have second 
order weak derivatives, they are very hard (if all possible) to identify. On the other 
hand, since strong moment solutions do have second order weak derivatives, naturally 
they are expected to satisfy the PDE \1.2\ almost everywhere in and to fulfill the 
boundary condition (2.3) pointwise on dft (cf fS^ [63 f ). 



3. Discretization and solution methods. The vanishing moment method 



reduces the problem of solving (1.2), (2. 3) to a problem of solving (2.2),(2.4)i,(2.5)i 
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for each fixed e > 0. Since (2.2) is a nonhnear biharmonic equation, one can use any of 
wealthy amount of existing numerical methods for biharmonic problems to discretize 
the equation. Although other types of numerical methods are applicable, here we 
focus on Galerkin type methods such as finite element methods, mixed finite element 
methods, discontinuous and spectral Galerkin methods [18l |9l |22j [8] . Throughout this 
section, we assume n = 2, 3. 

3.1. Finite element methods in 2-d. In the two-dimensional case many finite 
element methods, such as confirming Argyris, Bell, Bogner-Fox-Schmit and Hsieh- 
Clough-Tocher elements and nonconforming Adini, Morley, and Zienkiewicz elements, 
were extensively developed in 60 's and 70 's for the biharmonic problems. A beautiful 
theory of plate finite element methods was also established (cf. [8]). Naturally, one 



would want to solve problem (2.2),(2.4)i,(2.5)i by using and adapting these well- 
known plate finite element methods. That is exactly what we are going to do next. 
For the sake of presentation clarity, here we only discuss the confirming finite element 
methods, and refer to [63 for a detailed development of nonconfirming finite element 
methods for problem (l2l!2| , ([23 1 , (l2^i. 



The variational formulation for (2.2),(2.4)i, (|2.5|)i is defined as: Find G H^{Q) 



with = g a.e. on dQ such that for any v G H'^{Q) H Hq{Q) there holds 

-e{Au',Av) + {F{D^u',\/u',u',x),v) = -(^^ |^)- (3-1) 

Let 7^ be a quasiuniform triangular or rectangular mesh with mesh size h G (0, 1) for 
the domain Q C R^. Let Ug C H'^{Q) denote one of confirming finite element spaces 
(as mentioned above) whose functions take the boundary value g at all nodes on dCl. 
Then our finite element method is defined as: Find uf^ G Ug such that 

-s{Aul,Avh) + {F{D^ul,Vul,ul,x),Vh) = -(e^ G U^}. (3.2) 

In ^ we shall present several numerical experiment results for the Monge- Ampere 
type equations to show the excellent performance of the Argyis finite element method. 
Convergence and error analysis of the above scheme and other finite element schemes 
will be presented in forthcoming papers (also see [63 ). 

3.2. Mixed finite element methods in 2-d and 3-d. Along with the theory 
of plate finite element methods, another beautiful theory of mixed finite element 
methods was also extensively developed in '70s and '80s for the biharmonic problems 
in 2-d (cf. \9, TSl |35]). It is interesting to point out that all these 2-d mixed finite 
element methods can be easily generalized to solving 3-d biharmonic problems and 
general fourth order quasilinear PDEs (cf. [32l [4TJ [42] ) . 

Because the Hessian matrix D^u^ appears in ( |2.2[ ) in a nonlinear fashion, to design 
a mixed method we are "forced" to introduce cr^ := D^u^ (not := Au^ alone) 
as additional variables so the mixed method simultaneously seeks and <j^. This 
observation then excludes the usage of the popular family of Ciarlet-Raviart mixed 
finite element methods (originally designed for the biharmonic problems) [181 [HI ? 
the other hand, the observation suggests to try Hermann- Miyoshi mixed elements 
[35l [53| [6Ql |6T| ^8^ and Hermann- Johnson mixed elements [35j [531 ES] both use cr^ as 
additional variables. 

To define Hermann- Miyoshi type mixed finite element methods, we first de- 



rive the following mixed variational formulation for problem (2. 2), (2. 4 )i, (2. 5)3: Find 
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{u^.cr^) eVg xWs such that 



i=l 



£(divo■^V^;)+ (F(a^V^^^^^^a;),^') = (/,«) e Vq, 



(3.3) 
(3.4) 



where r^, z = 1, 2, • • • , (n — 1) denote the (n — 1) tangential directions at each point 
on dft^ ^ denotes the tangential derivative of g along r^, and 



Vg: = {ve H\ay, vlon = g}, Vo := {v e H\n); vlen 
Wo: = {fie /Xy = fiji, tin ■ n\9n = O}. 



0}, 



Let be a quasiuniform triangular or rectangular mesh if n = 2 and be a 
quasiuniform tetrahedronal or 3-d rectangular mesh if n = 3 for the domain Q. Let 
yh ^ be the Lagrange finite element space consisting of continuous piecewise 

polynomials of degree k{> 2) associated with the mesh Th. Let 



yh ._ yh ^ y^^ yh ._ yh ^ y^^ ^^h ._ ^yhy^Xn ^ ^h ._ ^yh^nXn ^ 



h ._ 



Based on the variational formulation ( |3.3[ )-( [3^ we define our (Hermann-Miyoshi 
type) mixed finite element methods as follows: Find {uf^^al) G x such that 



{(Jh^lih) + (V<,div/i/,) = ^{-^^^l^h^-^i) "^l^h e W^, 



i=l 



(3.5) 
(3.6) 



Similarly, we can define variants of the above scheme as those proposed in [68] 
as well as Hermann- Johnson type mixed methods. In ^ we shall present several 
numerical experiment results for the above scheme applying to the Monge- Ampere 
type equations. Convergence and error analysis of the above scheme and other mixed 
finite element schemes will be presented in a forthcoming paper (also see [63]). 

Remark 3.1. Besides the finite element and mixed finite element discretiza- 
tion methods, one can also approximate problem (2.2)^(2.4)i^(2.5)3 by discontinuous 
Galerkin methods '2l [H 



\3l\\M\SM (^'^d spectral Galerkin methods |3 [771[gP| /. 

It should be pointed out that these methods are dimens on-independent, hence, can be 
used in both 2-d and 3-d cases. We refer to 163] for a detailed exposition. 

3.3. Remarks on second order fully nonlinear parabolic equations. By 

adopting the method of line approach, generalizations of the numerical methods dis- 
cussed in previous subsections to the corresponding parabolic equations (4.11) and 
(4.12) are standard (cf. [32l [38] and references therein). Assuming that an implicit 
time stepping method such as the backward Euler and the Crank-Nicolson schemes 
will be used for time discretization, then at each time step we only need to solve a 



fully nonlinear elliptic equation of the form (2.2). As a result, all numerical methods 



discussed in §3 . lf^3^ immediately apply. On the other hand, it should be pointed out 
that the convergence and error analysis of all fully discrete schemes are expected to be 
harder, in particular, establishing error estimates which depend on polynomially 
instead of exponentially will be very challenging (cf. [40] [H] |42l |39l |44] |43]). 
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3.4. Remarks on nonlinear solvers and preconditioning. After equations 



(2.2) and (4.11) are discretized by any of above discretization methods, we get a 
strong nonlinear algebraic system to solve. To the end, one has to use one or another 
iterative solution method to do the job. In all numerical experiments to be given in 
^ we use preconditioned Newton iterative methods as our nonlinear solvers. A few 
fixed point iterations might be needed to generate initial guess for Newton type iter- 
ative methods. Another strategy which we are currently investigating is the following 
"multi-resolution" strategy: first compute a numerical solution using a relatively large 
£, then use the computed solution as an initial guess for the Newton method at the 
finer resolution e. Regarding to preconditioning, we use the simple ILU preconditioner 
in all simulations of ^ We plan to use more sophisticate multigrid and Schwarz (or 
domain decomposition) preconditioners when the amount of computations becomes 
intensive and large in 3-d. 

4. Applications. In this section, we shall apply the vanishing moment method- 
ology outlined in the previous section to several classes of specific second order fully 
nonlinear PDEs. 

4.1. Monge- Ampere type equations. Monge- Ampere type equations refer to 
a class of second order fully nonlinear PDEs of the form (cf. (131 Ull HZl US [51]) 

F{D''u\ Du\u\ x) := det{D^u'') - /(V^l^ x) = 0, (4.1) 

Note that from now on we shall always use u'^ to denote a solution of a second order 



fully nonlinear PDE we intend to solve. Equation (4.1 ) reduces to the classical Monge- 
Ampere equation 

det{D^u^) = f{x) 

if fiVvP^u^^x) = f{x) > 0, and to Gauss curvature equation 

det{D^u'')=K{l^\\/u''\^)^ 

if f{Vu^,u^,x) = K{1 + iV^i^p) 2 . Where the constant K is a prescribed Gauss 
curvature. Monge- Ampere type equations are the best known second order fully 
nonlinear PDEs, they arise in differential geometry and applications such as mass 
transportation and meteorology. It is well-known that Monge- Ampere type equations 
are elliptic only in the set of convex functions (cf. [Ml [51]). So their viscosity solutions 
are defined as convex functions in the sense of Definition 12.21 



The vanishing moment approximation (2.2) to (4.1) reads as: 

- eA^u' + det{D^u') = f{\/u', u',x) (e > 0). (4.2) 
For each fixed s > 0, this is a quasilinear fourth order PDE with Hessian type non- 



linearity. It is complemented by boundary conditions (2.4)i and (|2.5[)i (or (2.5)2). 



For the classical Monge- Ampere equation, it can be shown that [36l [63] 



• For each fixed £ > 0, problem ( [42| ), ([2^i, ( [23| )i (or ([2^2) has a unique 
weak solution in W^''^{Q) for n = 2,3. 

• det{D'^u^) > and Au^ > in 1] for £ > when n = 2, 3. 

• is convex in 1] for e > when n = 2. 

• The Dirichlet problem for the classical Monge- Ampere equation has a unique 
convex strong moment solution, which coincides with the unique convex vis- 
cosity solution of the same problem when n = 2. 
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The above results immediately imply that in the two dimensions the vanishing moment 
method indeed works for the classical Monge- Ampere equation and the notion of 
moment solutions and the notion of viscosity solutions are equivalent in this case. 



Remark 4.1. Recall that when n = 2, the Dirichlet problem (1.4) -(1.5) has at 

most two solutions (see 1231). An amazing numerical discovery to he given in ^^is 

that if we restrict e in (4.2) to e then lim also exists and the limit is nothing 
^0- 



hut the other solution solution of problem (1.4)-(1.5) which is concave (see Figures 



5.1-5.6)! 



4.2. Pucci's equations. Pucci's extremal equations are referred to the following 
two families of fully nonlinear PDEs (cf. j46l [T3]) 



Ma,[u] := aAu + (1 - na)\n{D'^u^) = f{x), 
ma[u] := aAu + (1 - na)Xi{D'^u^) = f{x) 



(4.3) 
(4.4) 



for < a < ^. Where Xn{D'^u^) and Xi{D'^u^) denote the maximum and minimum 
eigenvalues of the Hessian matrix D'^u^ . In the 2-d case, the above equations can be 
rewritten in terms of AvP and dei^D'^vP) (cf. [31 ). 

The vanishing moment approximations to (|4.3[) and (|4.4|) are defined as 



■ sA^u' 



■aAu^ 
■aAu^ 



(l-na)Xn{D^u') = f{x), 
(l-na)Xi{D'^u^) = f{x), 



(4.5) 
(4.6) 



which should be complemented by boundary conditions (|2.4[)i and (2.5)i (or (|2.5[)2)» 



4.3. Infinite Laplace equation. The infinite Laplace equation refers to the 
following degenerate quasilinear PDE: 



F{D^u^, Du^.u^, x) := A^oU^ = 0, 



(4.7) 



where 







AoqvP can be regarded as the limit of the p-Laplacian ApU^ := div(|Vii^|^ '^Vu^) as 
p ^ oo, it also can be derived as the Euler-Lagrange equation of the functional 



I{v) := ess sup^^^|Vv(x)|, 



whose minimizers are often called "absolute minimizers" [3 . Besides its mathematical 
appeals, the infinite Laplace equation also arises from image processing, geography, 
and geology applications [4, 15 . Although the infinite Laplace equation is only a 
degenerate quasilinear PDE, not a fully nonlinear PDE, it is very difficult to solve 
numerically. This is because the infinite Laplace equation does not have classical 
solutions in general [3 , and since it is not in divergence form, its weak solutions are 
defined and understood in the viscosity sense. We refer to [64 for recent developments 
on finite difference approximations of the infinite Laplace equation. 



Here we propose the following vanishing moment approximation for (4.7): 



(4.8) 



which is complemented by boundary conditions (2.4)i and (2.5)i (or (2.5)2). In ^ we 
shall present numerical results which show that the vanishing moment approximation 
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exactly converges to the unique viscosity solution of the Dirichlet problem for (4.7). 
This is another example which shows that the notion of moment solutions and the 
notion of viscosity solutions coincide. 

Remark 4.2. It is easy to see that the above vanishing moment method also 
applies to the p-Laplacian equation —ApU^ = f for 1 < p < oo. 

4.4. Second order fully nonlinear parabolic PDEs. We first like to note 
that there are several different versions of legitimate parabolic generalizations to ellip- 



tic PDE (1.2) (cf. ^l[70 ). In this paper, we shall only consider the following widely 
studied (and it turns out to be the "easiest") class of second order fully nonlinear 
parabolic PDEs: 

F{D^u^, V^x^ u^, X, t) - u^^ = 0, (4.9) 
assuming that F{D'^u^ ^VvP ^v!^ ^x^t) is elliptic. Clearly, this is the most natural 



parabolic generalization to equation (1.2). For example, the corresponding parabolic 
Monge- Ampere type equation reads as 

dei{D^u^) -u^t= /(V^i^ x, t) > 0. (4.10) 

In the past two decades the viscosity solution theory has been well developed for equa- 



tions (4.9) and (4.10), see [571 [TQl |52] . On the other hand, numerical approximation 
to these fully nonlinear parabolic PDEs is a completely untouched area. To the best 
our knowledge, no numerical result (in fact, no attempt) is known in the literature. 
Similarly, we can define the vanishing moment method and the notion of moment 



solutions for initial and initial-boundary value problems for (4.9), and then ask the 
same questions as we did in §2.2| We leave this as an exercise to interested readers 
and refer to [63 for a detailed exposition. 

Following the derivation of §2.2[ we propose the following vanishing moment ap- 



proximations to (|4.9[) and (4.10), respectively. 



F{D^u', V^i^ x, t) - sA^u' - = 0, (4.11) 
det{D^u') - sA^u' -ul = f(yu', ^i^ x, t), (4.12) 

each of the above equations is a fourth order quasilinear parabolic PDEs. 

5. Numerical experiments. In this section, we shall present a number of nu- 
merical experiment results obtained by using the vanishing moment method together 
with the numerical methods proposed in ^ Both 2-d and 3-d tests will be pre- 
sented. All the 3-d tests are obtained by a Hermann-Miyoshi type mixed finite el- 
ement method, while the 2-d tests are computed by using both the Argyris (plate) 
finite element method and the Hermann-Miyoshi mixed finite element method. 

5.1. Two-dimensional numerical experiments. The numerical solutions of 
the first seven tests are computed using the Argyris finite element method. 



Test 1: In this test we solve the Monge- Ampere problem (1.4)-(1.5) on the unit 
square Vt = (0, 1)^ with the following data: 

fix.y) 1, gix.y) 0. 

We remark that problem ([L4|-([L5| has a unique convex viscosity solution but does 
not have a classical solution (cf. I31j). 
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(2.4) 



Recall that the vanishing moment approximation of (1.4)-(1.5) is problem (4.2), 
(2.5 )i with the above / and g. We discretize problem (4.2), (2.4)i, (2.5 )i 
using the Argyris plate element as described in §3.1| Figure \5A\ displays the computed 
(moment) solutions using e = 10~^ (left graph) and s = — 10~^ (right graph). Clearly, 
the vanishing moment approximations correctly capture the convex viscosity solution 
(left graph) and the concave viscosity solution (right graph). Hence, the moment 
solutions coincide with the viscosity solutions (see [36 for a rigorous proof). 





i 



i 



Fig. 5.1. Computed (moment) solutions of Test 1: Graph on left corresponds to e = 10 



and graph on right corresponds to e ■ 



-10- 



To have a better view of the convexity of the computed solution, we also plot 
selected cross sections of the left figure in Figure [5?T| The cross sections clearly show 
that the computed solution is a convex function. In particular, there is no visible 
smear at the boundary. 
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0.5 0,6 0.7 



Fig. 5.2. x-cross sections (left figure) of the left graph in Fig. 
0.1,0.3,0.5,0.7,0.9 (indicated respectively by asterisk, circl e, pl us sign, 
gle); y-cross sections (right figure) of the left graph in Fig. 



5.1 



5 A 



(indicated respectively by asterisk, circle, plus sign, squre, and triangle) 



at X — 
squre, "and trian- 
aty = 0.1,0.3,0.5,0.7,0.9 



Test 2: The only difference between this test and Test 1 is that the datum 
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functions are now chosen as 



fix, y) = il + (x^ + y^))e^^'+y'\ g{x, y) = { 



e.V2 


if 


X 


= 0, 




if 


y 


= 0, 


e(i+^')/2 


if 


y 


= 1, 


e(i+!/')/2 


if 


X 


= 1, 



SO that vP{x^y) = ^e^^ +^ ^ is an exact solution of problem (1.4)-(1.5). Clearly, is 



a convex function, hence must be the unique convex viscosity solution of problem 



(|L4|)-(|L5| (cf^^ [51J. 



Figure 
and £ = — 10 



5.3 



shows the computed (moment) solutions using e = 10-^ (left graph) 
(right graph). Again, the vanishing moment approximations correctly 
capture the convex viscosity solution (left graph) and the concave viscosity solution 
(right graph), hence the moment solutions coincide with the viscosity solutions (see 
[36 for a rigorous proof). 





1 



- -1.1 



Fig. 5.3. Computed (moment) solutions of Test 2: Graph on left corresponds to e 
and graph on right corresponds to e = — 10~^. 



10" 



Test 3: Similar to Test 2, the only difference between this test and Test 1 is that 
the datum functions are now chosen as 



so that u^{x^ y) 





if X 


= 0, 




if y 


= 0, 


^(l + x2)t 


if y 


= 1, 


2^(l + y2)t 


if X 


= 1, 



2V2 

3 



9{x,y) 



(x^ + y^)^ is the unique convex viscosity solution of problem 



Figure 
and e 



5.4 



10- 



displays the computed (moment) solutions using e = 10~^ (left graph) 
^ (right graph). As expected, the vanishing moment approximations 
correctly capture the convex viscosity solution (left graph) and the concave vis- 
cosity solution (right graph), hence the moment solutions coincide with the viscosity 
solutions (see [36 for a rigorous proof). 

Test 4: Again, the only difference between this test and Test 1 is that the datum 
functions are now chosen as 

fi^^y) = (i-x-yf g = o. 
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On the other hand, mathematicahy there is a significant difference between these two 
test problems. Note that /(x, y) = on the line x -\- y = 1 in the domain Q = (0, 1)^. 
Hence, problem (1.4)-([L5]) is known as a degenerate Monge- Ampere problem (cf. 

m)- 



Figure 5.5 displays the computed (moment) solutions using e = 10~^ (left graph) 
and e = — 10~^ (right graph). Once again, the vanishing moment approximations 
correctly capture the convex viscosity solution (left graph) and the concave viscosity 
solution (right graph), hence the moment solutions coincide with the viscosity solu- 
tions (see [36 for a rigorous proof). In addition, our numerical result shows that the 
vanishing moment method is robust with respect to the degeneracy of the underlying 
PDE. 




Test 5: Once again, the only difference between this test and Test 1 is that the 
datum functions are now chosen as 

f{x,y) = -y^ g = 0. 



Mathematically, the difference between this test problem and Test 1 is even more 
dramatic because not only /(x, y) = on the line x — y = but also f changes sign 
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(hence the PDE changes type) in ft. To the best of our knowledge, there is no viscosity 
solution theory for this type Monge- Ampere problems in the literature. However, the 
vanishing moment method seem works well for this problem. Our numerical results 
indicate existence of both convex and concave moment solutions. 



Figure 5.6 displays the computed convex (moment) solution using e = 10 ^ (left 
graph) and the computed concave (moment) solution using e = — 10~^ (right graph). 



Surface: u Height: u 



Surface: u Height: u 





t 



Fig. 5.6. Computed (moment) solutions of Test 5: Graph on left corresponds to e — 10 
and graph on right corresponds to e — — 10~^. 



Again, to have a better view of the convexity of the computed solution, we also 



plot selected cross sections of the left figure in Figure |5.6| The cross sections clearly 
show that the computed solution is a convex function. In particular, there is no visible 
smear at the boundary. 
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5.6 at X 



Fig. 5.7. x-cross sections (left figure) of the left graph in Fig. 
0.1,0.3,0.5,0.7,0.9 (indicated respectively by asterisk, circl e, pl us sign, squre, and trian- 
gle); y-cross sections (right figure) of the left graph in Fig. 5.6 at y — 0.1,0.3,0.5,0.7,0.9 
(indicated respectively by asterisk, circle, plus sign, squre, and triangle) 



Test 6: In this test we solve the following Gauss curvature (or /C-surface) equation 
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(cf. [Ml [501 Ha]) 

det (D^u^) 



y 



in Q := (-0.57, 0.57)^ 
on dQ, 



(5.1) 
(5.2) 



where K > is a given constant Gauss curvature. Note that the above problem is a 
special case of problem (|4l]),(|23| with f{\/u^,u^,x,y) = K{1 + IV^/^P)"^, n = 2, 



and g{x, y) = x -\- y — 1. 

It was proved by Guan [48^ that there exists > such that for each K G 
[0,i^*) problem (5.1 )-( [5^ (with more general Dirichlet data) has a unique convex 
viscosity solution. Theoretically, it is very difficult to give an accurate estimate for 
the curvature upper bound K* . However, hand, this offers an ideal opportunity for 
numerical analysts to help and to contribute. It turns out that the vanishing moment 
method proposed in this paper works very well for such a problem, hence it might 
provide a useful tool and answer to the challenge. 

Since we are only interested in convex solutions of the Gauss curvature equation, 

displays the computed convex (moment) 



so we restrict e > in (4.2). Figure 5 

and K 



solution using s = 10 



0.1,1,2,2.1, respectively. We note that our com- 



Surface: u Height: 




Surface: u Height: u 





i 



Fig. 5.8. Computed (moment) solutions of Test 6: £ = 10 ^ and K — 0.1,1,2,2.1. 
Graphs are arranged row-wise. 



puter code stops producing a convergent numerical solution for K = 2.2. Hence we 
conjecture that ^ 2.1 for the above test problem. 



Test 7: In this test, we solve problem (4.7), (2.3) over the domain Q = (— ^, ^)^ 
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Fig. 5.9. Computed moment solution (left graph) and the exact viscosity solution (right 
graph) of Test 7. s = 10"^ 



with the fohowing boundary datum function 



4 

X3 — 



if X 
if y 
if X 
if y 



so that vPiyX^y) = x^ — y^ is the unique viscosity solution (cf. [3 ). 

We remark that this is an important example in the theory of absolutely minimiz- 
ing functions since is the least regular absolutely minimizing function known in the 
case of the Euclidean norm (cf. [3^ and the references therein). It is easy to check that 
is a Holder continuous function with exponent |. However, it is not twice differen- 
tiable on the axes. We also note that (see ^4.3|) the infinite Laplace equation (4.7) is 



only a second order (degenerate) quasilinear (instead fully nonlinear) PDE, however, 
the complicate and nondivergence structure makes the infinite Laplace equation very 
difficult to analyze theoretically and to compute numerically. 

Figure 5.9 displays the computed (moment) solution using e = 10-^ (left graph) 
and the exact solution (right graph). Once again, the vanishing moment ap- 
proximation correctly captures the viscosity solution u^, hence, the moment solution 
coincides with the viscosity solution (see [36 for a rigorous proof). 

The numerical solutions of the next two tests are obtained by using the Hermann- 
Miyoshi mixed finite element method with piecewise quadratic shape functions. 

Test 8: This test is a re-run of Test 1 but using the quadratic Hermann-Miyoshi 
mixed finite element method. Figure [5?Tq| is the counterpart of Figure [5?T| We remark 
that the mixed method also produces an approximation to the Hessian matrix D^u^ ^ 
which is not shown here. Clearly, the numerical results of Test 1 and Test 8 have the 
same accuracy. However, it should be noted that the mixed method runs about 20 
times faster than the Argyris method on this test problem. 

Test 9: This test solves, using the quadratic Hermann-Miyoshi mixed finite 
method, the Monge- Ampere problem (1.4)-(1.5) on the unit square Q = (0, 1)^ with 
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Max: 

I 



Fig. 5.10. Computed (moment) solutions of Test 8: Graph on left is the computed 



with £ = 10 ^ and graph on right is the computed u% with , 



-10" 





FiG. 5.11. Computed (moment) solutions of Test 9 by the quadratic Hermann- Miyoshi 
mixed method, e — . 



the following data: 



f{x,y) 



(4- 



so that = ^ — — is an exact (convex) solution. We note that problem 
( |1.4| )-( pr5| ) has e xact two solutions, one is convex and the other is concave (cf.[23 ). 

Figure 5.11 displays the computed (moment) solution using e = 10~^ (left 
graph) and its error (right graph). As expected, the vanishing moment approxima- 
tion correctly captures the convex viscosity solution u^. Hence the moment solution 
coincides with the viscosity solution (see [36 for a rigorous proof). Again, we remark 
that the mixed method also gives an approximation to the Hessian matrix D^u^, 
which is not shown here, and the mixed method runs about 20 times faster than the 
Argyris method for solving this test problem. 

5.2. Three-dimensional numerical experiments. In this subsection we present 
two numerical tests on computing moment (and viscosity) solutions of the Monge- 
Ampere problem (1.4)-(1.5) in the unit cube Q = (0, 1)^. Numerical approximations 



of fully nonlinear PDEs in 3-d is known to be very difficult. To the best of our knowl- 
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edge, no 3-d numerical results are given in the literature for the Monge- Ampere type 
fully nonlinear PDEs. 

Test 10: Consider the Monge-Ampere problem (1.4)-(1.5) on the unit cube 
(0, 1)^ with the following data: 



f{x, y, = (1 + + + x'^) exp i^- 



y 



9{x,y,z) = exp 



y 



-)■ 



It is easy to verify that = exp ( ^ ) is a unique exact (convex) solution. We 

compute this solution using the vanishing moment method combined with a general- 
ized Hermann- Miyoshi type mixed finite element method using linear shape functions. 
Figure [5?T2| displays color plots of five x-slices of the computed (moment) solution 



(left graph) and its corresponding error (right graph). Figure 5.13 displays color 



plots of five z-slices of the computed (moment) solution uf^ (left graph) its correspond- 
ing error (right graph). As expected, the vanishing moment approximation correctly 
captures the convex viscosity solution u^. Hence the moment solution coincides with 
the viscosity solution. 





Fig. 5.12. x-slices of the computed (moment) solution of Test 10 by a generalized linear 
Hermann- Miyoshi mixed method, e — 10~^. 
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Test 11: Our last numerical test solves the 3-dimensional generalization of the 
test problem in Test 1. That is, we assume u satisfies the Monge- Ampere problem 



([L4])-([L5]) mQ = (0, 1)^ with the data 

f{x,y,z) 1, g{x,y,z) 0. 

We remark that the above problem has a unique convex viscosity solution but does 
not have a classical solution (cf. [5ll|3T]). There is no explicit solution formula for 
the boundary value problem. 



Figure 5.14| displays color plots of x-slices (left graph) and z-slices (right graph) 



of the computed (moment) solution uf^ using a generalized linear Hermann-Miyoshi 
mixed method. Once again, the vanishing moment approximation correctly captures 
the convex viscosity solution u^. Hence the moment solution coincides with the vis- 
cosity solution. 




6. Conclusions. In this paper we introduce a new notion of weak solutions, 
called moment solutions^ through a constructive limiting process, called the vanish- 
ing moment method^ for second order fully nonlinear PDEs. The notion of moment 
solutions and the vanishing moment method are exactly in the same spirit as the orig- 
inal notion of viscosity solutions and the vanishing viscosity method proposed by M. 
Crandall and P. L. Lions in [26 for the Hamilton- Jacobi equations, which is based on 
the idea of approximating a fully nonlinear PDE by a higher order quasilinear PDE. 
We first present a general framework of the vanishing moment method and the notion 
of moment solutions in Q We then apply the general framework to several classes 
of PDEs including the Monge-Ampere type equations, Pucci's extremal equations, 
the infinite Laplace equation, and second order fully nonlinear parabolic PDEs. We 
then propose two classes of numerical methods to discretize the fourth order "regu- 
larized/perturbed" vanishing moment approximation equations. Finally, we present 
a number of numeral experiments using the vanishing moment methodology together 
with the proposed numerical methods to demonstrate convergence and effectiveness 
of the vanishing moment method, as well as the relationship between the notion of 
moment solutions and the notion of viscosity solution for second order fully nonlinear 
PDEs. 

This paper provides a practical and systematic methodology/approach, which can 
be backed by rigorous PDE and numerical theories, for approximating second order 
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fully nonlinear PDEs. As a by-product, the moment solution theory will provide some 
insights to our understanding of the viscosity solution theory, and might provide 
a logical and natural generalization/extension for the notion of viscosity solution, 
especially, in the cases where there are no theories or the existing viscosity solution 
theory fails (such as Monge- Ampere equations of sub elliptic and hyperbolic types 
[16], and systems of second order fully nonlinear PDEs.) 

Acknowledgment: The first author would like to thank his former colleague 
Professor Bo Guan of Ohio State University for many stimulating discussions and 
helpful suggestions. 
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